##################################################################
# /* Author: Gautam Nair and Nicholas Sambanis */
# /* Violence Exposure and Ethnic Identification: Evidence from Kashmir */
# /* Figure 1: Self-Reported Violence Exposure and Ethnic Identification  */

##################################################################

# change working directory here
# setwd("")

##################################################################
# loading packages
##################################################################
library(Hmisc)
library(foreign)
library(sandwich)
library(lmtest)
library(numDeriv)
library(stargazer)
library(ggplot2)
library(plyr)
library(gridExtra)
library(dplyr)
library(scales)

##################################################################

rm(list=ls())

library(stargazer)

data.working <- read.csv("nei_kashmir_replication_dataset.csv")

dep.var <- c("identity_rank_indian_rev", 
"identity_choice_indian_both", 
"kashmir_status_india",  
"ind_pak_prefer_ind",
"bonus_donation_india", 
"ic_index_id_no_ind_pak",
"ic_index_id_ind_pak",
"protest_peaceful_endorse",
"protest_violent_endorse"
)

ind.violence <- c("violence_exp_index_dummy",
				"group2_violence",
				"group3_growth",
				"group4_int_inst",
                  "group5_geography")

ind.correlates <- c("age_18_27",               
                    "age_28_49",               
                    "gender_female",
                    "education_class12",
                    "religion_sunni",
                    "religion_shia",
                    "main_earner",
                    "index_income_wealth1",
                    "occupation_govt")

ind.districts <- c("constituency_1",  
                   "constituency_2", 
                   "constituency_3", 
                   "constituency_4", 
                   "constituency_5", 
                   "constituency_6", 
                   "constituency_7", 
                   "constituency_8", 
                   "constituency_9", 
                   "constituency_10", 
                   "constituency_11", 
                   "constituency_12", 
                   "constituency_13", 
                   "constituency_14", 
                   "constituency_15", 
                   "constituency_16", 
                   "constituency_17", 
                   "constituency_18", 
                   "constituency_19", 
                   "constituency_20", 
                   "constituency_21", 
                   "constituency_22")
# constituency X is excluded

ind.1 <- ind.violence
ind.2 <- c(ind.violence, ind.correlates)
ind.3 <- c(ind.violence, ind.correlates, ind.districts)

reg.models <-vector("list", length(dep.var)) 
se.models <-vector("list", length(dep.var))
matrix <- vector("list", length(dep.var))
x=0

for(i in 1:length(dep.var)){
  
  my.formula <-paste(dep.var[i],'~', paste(ind.1, collapse= ' + '))
  temp.model <- lm(my.formula, data=data.working)
  temp.se <- coeftest(temp.model, vcov=vcovHC(temp.model,type="HC2"))[,2]
  x = x+1
  se.models[[x]] <- temp.se
  reg.models[[x]] <- temp.model
  
  my.formula <-paste(dep.var[i],'~', paste(ind.2, collapse= ' + '))
  temp.model <- lm(my.formula, data=data.working)
  temp.se <- coeftest(temp.model, vcov=vcovHC(temp.model,type="HC2"))[,2]
  x = x+1
  se.models[[x]] <- temp.se
  reg.models[[x]] <- temp.model	
  
  my.formula <-paste(dep.var[i],'~', paste(ind.3, collapse= ' + '))
  temp.model <- lm(my.formula, data=data.working)
  temp.se <- coeftest(temp.model, vcov=vcovHC(temp.model,type="HC2"))[,2]
  x = x+1
  se.models[[x]] <- temp.se
  reg.models[[x]] <- temp.model	
  
}


##################################################################
# Chart
##################################################################

coef.list<- as.numeric(rbind(
summary(reg.models[[1]])$coefficients[2,1],
summary(reg.models[[2]])$coefficients[2,1],
summary(reg.models[[3]])$coefficients[2,1],
summary(reg.models[[4]])$coefficients[2,1],
summary(reg.models[[5]])$coefficients[2,1],
summary(reg.models[[6]])$coefficients[2,1],
summary(reg.models[[7]])$coefficients[2,1],
summary(reg.models[[8]])$coefficients[2,1],
summary(reg.models[[9]])$coefficients[2,1],
summary(reg.models[[10]])$coefficients[2,1],
summary(reg.models[[11]])$coefficients[2,1],
summary(reg.models[[12]])$coefficients[2,1],
summary(reg.models[[13]])$coefficients[2,1],
summary(reg.models[[14]])$coefficients[2,1],
summary(reg.models[[15]])$coefficients[2,1],
summary(reg.models[[16]])$coefficients[2,1],
summary(reg.models[[17]])$coefficients[2,1],
summary(reg.models[[18]])$coefficients[2,1],
summary(reg.models[[19]])$coefficients[2,1],
summary(reg.models[[20]])$coefficients[2,1],
summary(reg.models[[21]])$coefficients[2,1],
summary(reg.models[[22]])$coefficients[2,1],
summary(reg.models[[23]])$coefficients[2,1],
summary(reg.models[[24]])$coefficients[2,1],
summary(reg.models[[25]])$coefficients[2,1],
summary(reg.models[[26]])$coefficients[2,1],
summary(reg.models[[27]])$coefficients[2,1]
))

se.list<-
as.numeric(rbind(
se.models[[1]][2],
se.models[[2]][2],
se.models[[3]][2],
se.models[[4]][2],
se.models[[5]][2],
se.models[[6]][2],
se.models[[7]][2],
se.models[[8]][2],
se.models[[9]][2],
se.models[[10]][2],
se.models[[11]][2],
se.models[[12]][2],
se.models[[13]][2],
se.models[[14]][2],
se.models[[15]][2],
se.models[[16]][2],
se.models[[17]][2],
se.models[[18]][2],
se.models[[19]][2],
se.models[[20]][2],
se.models[[21]][2],
se.models[[22]][2],
se.models[[23]][2],
se.models[[24]][2],
se.models[[25]][2],
se.models[[26]][2],
se.models[[27]][2]
))

model.name <- rep(c("1 No Controls","2 Individual-level Covariates", "3 Individual-level Covariates + Constituency Fixed Effects"), 9)
model.number <- rep(c(1,2,3), 9)
variable.name <- c(
rep("1 Ranked Indian Identity (Reversed 1 to 4 )", 3), 
rep("2 Indian or Indian+ Kashmiri/Kashmiri Only (1/0)", 3),
rep("3 Kashmir Status: Remain Part of India (1/0)", 3),
rep("4 India/Pakistan: Prefer India (1/0)", 3),
rep("5 Donation to All India NGO (%)",3),
rep("6 Identity Index I", 3),
rep("7 Identity Index II", 3),
rep("8 Peaceful Protest Participation (1 to 5)", 3),
rep("9 Violent Protest Endorsement (1 to 5)", 3)
)  

# plotdata <- as.data.frame(cbind(coef.list,se.list,model.name,variable.name) )
plotdata <- as.data.frame(cbind(coef.list,se.list))
names(plotdata)[1] <- "coef.list"
names(plotdata)[2] <- "se.list"

plotdata$coef.list <- as.numeric(plotdata$coef.list)
plotdata$se.list <- as.numeric(plotdata$se.list)

plotdata
interval1 <- -qnorm((1-0.9)/2)
interval2 <- -qnorm((1-0.95)/2)

plotdata$lower.90 <- plotdata$coef.list-plotdata$se.list*interval1
plotdata$upper.90 <- plotdata$coef.list+plotdata$se.list*interval1
plotdata$lower.95 <- plotdata$coef.list-plotdata$se.list*interval2
plotdata$upper.95 <- plotdata$coef.list+plotdata$se.list*interval2

plotdata$ymin <- -(plotdata$se.list*2.5)
plotdata$ymax <- +(plotdata$se.list*2.5)
plotdata$model.name <- factor(model.name, levels = c("1 No Controls", "2 Individual-level Covariates", "3 Individual-level Covariates + Constituency Fixed Effects"))

plotdata$variable.name <- as.factor(variable.name)

plotdata$depvar <- as.factor("depvar")

plot.list <-vector("list", length(dep.var)) 
temp.font <- 9
temp.legend <-9
z=0
for(i in c(1,4,7,10,13,16,19, 22, 25)){
	q <- i+2
	tempdata <- plotdata[i:q,]
	if(i!=22){
		y.min <- -tempdata$se.list[1]*7
		y.max <- +tempdata$se.list[1]*7
	}
	if(i==22){
		y.min <- -tempdata$se.list[1]*10
		y.max <- +tempdata$se.list[1]*10
	}
	plot1 <- ggplot(tempdata, aes(shape=model.name, colour=model.name))
	plot1 <- plot1 + geom_linerange(aes(x = depvar, ymin = lower.90, ymax = upper.90), lwd = 0.75, position = position_dodge(width = 1/2))                                                     
	plot1 <- plot1 + geom_pointrange(aes(x = depvar, y= coef.list, ymin = lower.95, ymax = upper.95), lwd = 0.50, position = position_dodge(width = 1/2))
	plot1 <- plot1 + coord_flip()                             
	plot1 <- plot1 + theme_bw((base_size = temp.font))
	plot1 <- plot1 + geom_hline(yintercept = 0)
	plot1 <- plot1 + theme(axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank(),
							legend.title=element_blank()) + guides(fill = guide_legend(reverse=TRUE))
	plot1 <- plot1 + facet_wrap(~variable.name)	
	plot1 <- plot1 + scale_y_continuous(limits = c(y.min, y.max))
	z <- z+1
	plot.list[[z]] <- plot1
}


library(ggplot2)
library(gridExtra)
library(grid)


grid_arrange_shared_legend <- function(..., ncol = length(list(...)), nrow = 1, position = c("bottom", "right")) {

  plots <- list(...)
  position <- match.arg(position)
  g <- ggplotGrob(plots[[1]] + theme(legend.position = position))$grobs
  legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
  lheight <- sum(legend$height)
  lwidth <- sum(legend$width)
  gl <- lapply(plots, function(x) x + theme(legend.position="none"))
  gl <- c(gl, ncol = ncol, nrow = nrow)

  combined <- switch(position,
                     "bottom" = arrangeGrob(do.call(arrangeGrob, gl),
                                            legend,
                                            ncol = 1,
                                            heights = unit.c(unit(1, "npc") - lheight, lheight)),
                     "right" = arrangeGrob(do.call(arrangeGrob, gl),
                                           legend,
                                           ncol = 2,
                                           widths = unit.c(unit(1, "npc") - lwidth, lwidth)))
  grid.newpage()
  grid.draw(combined)

}

grid_arrange_shared_legend(plot.list[[1]], plot.list[[2]], plot.list[[3]], plot.list[[4]], plot.list[[5]], plot.list[[6]], plot.list[[7]], plot.list[[8]], plot.list[[9]],ncol = 1, nrow = 9)

png("tf_f_01_violence_id_correlation.png", height = 1800, width = 1400, units="px")

grid_arrange_shared_legend(plot.list[[1]], plot.list[[2]], plot.list[[3]], plot.list[[4]], plot.list[[5]], plot.list[[6]], plot.list[[7]], plot.list[[8]], plot.list[[9]],ncol = 1, nrow = 9)

dev.off()







